!************************************************************
! this subrouitne calculate point immediately downstream of
! shock in isentropic supersonic flow
!************************************************************
! Created by : B. G.
! Date       : 2015-08-13
! Revised    :
!************************************************************
! Nomenclature: 
!   ls   slope of shock line 14
!   p4s  p4 calculated from shock relationships
!   p4lrc p4 calculated from l-c-line(24) equations
!   dp   different of p4s and p4lrs, p4lrc-p4s
!   Vp~yp properties of functions of l-c-line
!   lp~Tp coefficients of functions of l-c-line
!   temp1 (g+1)/(2*g)
!   temp2 (g-1)/(2*g)
!   temp3 slope of secant method, (dp-dpf)/(eps4-eps4f), Eq. (g) of 17.4(g)
!   atemp temporary speed of sound
!   loc3  coefficient for locate point 3
!   dpf   dp of previous loop
!   eps4f eps4 of previous loop
!   i     loop indicator
!   
!************************************************************
! Warning :
!   1. angle in rad, not degree
!************************************************************

subroutine ShockP()
   use VariableDef
   implicit none
   ! given variables
   ! real*8::Vuu,pu,rhou,Tu,Mu,x1,y1,V1,theta1,p1,rho1,eps1,x2,y2,V2,theta2,p2,rho2
   ! temporary variables
   real*8::ls,p4s,p4lrc,dp
   real*8::Vp,pp,rhop,thetap,Mp,alphap,yp
   real*8::lp,Qp,Sp,Tp
   real*8::temp1,temp2,temp3,atemp,loc3
   real*8::dpf,eps4f
   integer::i
   ! calculate variables
   ! real*8::x4,y4,V4,p4,rho4,theta4,Vu4,Vv4,beta
   
   eps4=eps1
   temp1=(g+1)/(2*g)
   temp2=(g-1)/(2*g)
   i=1
   do
      ! determine p4s, theta4, V4 and rho4
      call Shock()
      p4s=p4
      ! determine ls
      ls=tan((eps1+eps4)/2)

      Vp=(V2+V4)/2
      pp=(p2+p4)/2
      rhop=(rho2+rho4)/2
      thetap=(theta2+theta4)/2
      call Thermo(pp,rhop,Vp,atemp,T4,Mp)
      alphap=asin(1/Mp)
      ! determine lp
      lp=tan(thetap+alphap)
      ! determine x4 and y4
      x4=(-ls*x1+lp*x2+y1-y2)/(lp-ls)
      y4=y1+ls*(x4-x1)
      yp=(y2+y4)/2
      ! determine coefficients along l-c-line
      Qp=sqrt(Mp**2-1)/(rhop*Vp**2)
      Sp=delta*sin(thetap)/(yp*Mp*cos(thetap+alphap))
      Tp=-Sp*(x4-x2)+Qp*p2+theta2
      p4lrc=(Tp-theta4)/Qp
      dp=p4lrc-p4s
      ! justify if point 4 arrived
      if (abs(dp) < ep) exit
      if (i==1) then
         eps4f=eps4
         dpf=dp
         eps4=asin(sqrt(temp1*p4lrc/pu+temp2)/Mu)   ! from Eq. 17.90
      else
         temp3=(dp-dpf)/(eps4-eps4f)   ! from Eq. (g) of 17.4(g)
         eps4f=eps4
         dpf=dp
         eps4=eps4-dp/temp3   ! from Eq. (g) of 17.4(g)
      end if
      i=i+1
   end do   
   call Thermo(p4,rho4,V4,atemp,T4,M4)
   Vv4=V4*sin(theta4)
   Vu4=V4*cos(theta4)
   ! locate point 3 midway between point 2 and 4
   loc3=0.5
   x3=x2+loc3*(x4-x2)
   y3=y2+loc3*(y4-y2)
   V3=V2+loc3*(V4-V2)
   p3=p2+loc3*(p4-p2)
   theta3=theta2+loc3*(theta4-theta2)
   rho3=rho2+loc3*(rho4-rho2)
   
end subroutine ShockP
